Kinetic theory for systems of self-propelled particles with metric-free interactions 
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A model of self-driven particles similar to the Vicsek model [Phys. Rev. Lett. 75 (1995) 1226] 
but with metric-free interactions is studied by means of a novel Enskog-type kinetic theory. In this 
model, N particles of constant speed Vo try to align their travel directions with the average direction 
of a fixed number of closest neighbors. At strong alignment a global flocking state forms. The 
alignment is defined by a stochastic rule, not by a Hamiltonian. The corresponding interactions are 
of genuine multi-body nature. The theory is based on a Master equation in 3N-dimensional phase 
space, which is made tractable by means of the molecular chaos approximation. The phase diagram 
for the transition to collective motion is calculated and compared to direct numerical simulations. 
A linear stability analysis of a homogeneous ordered state is performed using the kinetic but not 
the hydrodynamic equations in order to achieve high accuracy. In contrast to the regular metric 
Vicsek-model no instabilities occur. This confirms previous direct simulations that for Vicsek-like 
models with metric-free interactions, there is no formation of density bands and that the flocking 
transition is continuous. 
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I. INTRODUCTION 

One of the most important unsolved problems in sta- 
tistical physics is finding a global description of far-from- 
cquilibrium systems with many interacting objects. Such 
a unified theory would be especially useful for biological 
systems since they operate far from thermal equilibrium. 
Instead of looking for such a general theory we focus on a 
minimal nonequilibrium model which still displays inter- 
esting physics such as pattern formation and collective 
motion. The goal is to provide inspiration for a more 
general approach by constructing a quantitative theoret- 
ical framework for the minimal model. Wc consider a 
model similar to the Vicsek-model (VM) of self-propelled 
particles [l], Q which is simple enough to be treated nu- 
merically and analytically. The VM was introduced in 
1995 to describe the swarming of fish and birds 0]. In 
this model, pointlikc particles are driven with a constant 
speed. At each time step, a given particle assumes the 
average direction of motion of the particles in its neigh- 
borhood, with some added noise. This model constitutes 
a dynamical version of the 2D XY model, because the ve- 
locity of the "bird" , like the spin of the XY model, also 
has fixed magnitude and continuous rotational symme- 
try. As the amplitude of the noise decreases, the system 
undergoes a phase transition from a disordered state, in 
which the particles have no preferred global direction, 
to an ordered state, in which the particles move col- 
lectively in the same direction. Hence, unlike the XY 
model, Vicsek's model exhibits long-range orientational 
order at non-zero noise. This surprising fact motivated 
renormalization group studies by Toner and Tu [|[ which 
confirmed the stabilization of the ordered phase far from 
the flocking threshold. The phase transition was origi- 
nally thought to be continuous Q, but recent numeri- 



cal work [5( indicates that the transition is discontinuous 
with strong finite-size effects. The numerical studies also 
revealed that large density waves develop right next to 
the threshold while still maintaining global orientational 
order. 

Recently, it was shown by means of an Enskog-like ki- 
netic theory that the orderedphase of the VM is linearly 
unstable near the threshold [6J . This instability was pro- 
posed as a possible explanation for the density waves and 
discontinuous nature of the phase transition. Another 
study for a related model with continuous time found a 
similar instability by means of a Boltzmann equation Q , 
see also Ref. @. The VM assumes interaction with all 
neighbors within a fixed metric distance. Recent experi- 
ments by Ballcrini et al @ on flocks of several thousand 
starlings indicate that this interaction rule might not be 
appropriate for animal flocks. Instead, it was discovered 
that each bird interacts on average with a fixed num- 
ber of neighbors, typically six to seven. This constitutes 
a topological or metric-free interaction because not the 
metric distance is relevant but who are the closest neigh- 
bors. Ballerini et al argue further that due to evolution- 
ary pressure the main goal of interaction among individ- 
uals is to maintain cohesion. By comparing simulations 
with the regular VM and a modified VM with metric-free 
interactions they found that flocks, when facing preda- 
tors, kept cohesion much better in the metric-free model. 
This further supports the idea that metric-free interac- 
tions should be dominant in animal flocks. While quite a 
number of analytical and numerical studies about self- 
propelled particles with metric interactions have been 
published ToT- 18j , not much exists for topological inter- 



actions [L9 - 2L | . In particular, no rigorous theory for the 
metric-free model of Ref. @ exists. In order to construct 
a theory which can be applied directly to this minimal 
computer model as well as to experiments we adopted 
the original genuine multi-particle interactions and did 
not restrict ourselves to binary interactions. Since other 
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2D experiments on shoaling fish estimate the number of 
tracked neighbors to be between three to five [22| we ex- 
plored a range of interaction partner numbers between 
two and seven. 

The main results of this paper are a) the rigorous 
derivation of an Enskog-likc kinetic equation for the one- 
particle density for the metric-free model of Ref. [|[ from 
first principles, and b) a linear-stability analysis of this 
kinetic equation, which showed that the flocking state is 
linearly stable against perturbations of any wavelength. 

The remainder of this paper is structured as follows. 
In Sec. II we define the metric-free model. In Sec. Ill 
we set up an exact equation for the TV-particle proba- 
bility density and derive the kinetic equation Eq. (fT8| . 
In Sec. IV homogeneous solutions of this equation are 
discussed and the phase diagram of the order-disorder 
transition is calculated. Sec. V deals with the linear sta- 
bility analysis of the ordered phase, and Sec. VI describes 
direct simulations. A summary is given in Section VII. 
Details concerning hypergeometric functions, exact solu- 
tions for special cases and integral tables are relegated 
to Appendix A, B, and C, respectively. In Appendix D 
the Enskog kinetic approach is compared with the corre- 
sponding Boltzmann approximation. 



II. MODEL 

We consider a metric-free version of the VM, which 
was introduced in Ref. |9|. This two-dimensional model 
consists of N pointlikc particles with continuous spatial 
coordinates Vi(t) and velocities Vj(i) which evolve via two 
steps: streaming and collision. During a time step r, 
particles stream ballistically: x^(< + r) = Xi(t) + TVi(t). 
The magnitude of the particle velocities is fixed to vq. 
Only the directions Qi of the velocity vectors are updated 
in the collision step by first finding the M — 1 closest 
neighbors for a given particle i where M > 2 is a fixed 
parameter. The directions of motion of particle i and its 
neighbors determine the average direction, 
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The new flying directions follow as Oiit + r) = &i(t) 
where £j is a random number chosen with a uniform prob- 
ability from the interval [—r]/2, rj/2]. In the context of the 
VM P, Q this constitutes an angular noise model with 
noise strength rj. The particles are always updated in par- 
allel. Note, that in the original VM Q the number of in- 
teraction partners is fluctuating and density-dependent, 
whereas the range of interactions is fixed to the radius R 
of a circle around a given particle. It is the opposite in 
the metric-free model. Here, the number of interaction 
partners, M, is a fixed parameter but the interaction 
range fluctuates. For example, for small local density 
p(r), the next neighbors can be far away but no mat- 
ter how geographically isolated a particle is, it is always 



connected with M — 1 others via the metric-free interac- 
tion rule, Eq. dp To potentially allow comparison with 
experiments (9l. l23l| and simulations a large number 
M = 5 to 8 can be chosen. 

We define an effective interaction range i? c ff by inte- 
grating the density over a circle and equating the result 
with the partner number M: 
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resulting in 
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Another important length scale is the mean free path 
(mfp) given by the distance a particle travels between 
collisions, 



A = r v . 



(4) 



Note that the mfp is density-independent in VM-like 
models because of the discrete nature of the dynamics 
and because the particles have zero volume. 

The metric-free model is then characterized by four 
dimcnsionlcss control parameters: the noise strength 77, 
the ratio of the mfp to the effective interaction radius, 
A = X/R c g, the partner number M, and the normal- 
ized system size, L = L/X, of the L x L simulation box 
with periodic boundary conditions. Secondary parame- 
ters of less obvious physical relevance, such as the total 
particle number N and the average density po = N/L 2 
can be easily expressed in terms of the four character- 
istic parameters, for example, using Eq. ([3|) one finds, 
N = MA 2 L 2 /n. The thermodynamic limit, N — > 00, is 



then equivalent to L 
constant. 
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III. KINETIC THEORY 

Recently, a kinetic formalism for the Vicsek model be- 
yond the Boltzmann theory has been developed . Such 
an approach is particularly useful for the metric-free VM, 
where particles always interact with about four to seven 
neighboring particles at once. These genuine multi-body 
interactions cannot be described by the binary collision 
approximation of the Boltzmann equation plj . 

The starting point for the kinetic formalism is a 
discrete-time Master equation in 3N-dimensional phase 
space for the N-particle probability density 
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where XW = (xi,x 2 ,..., XA r) and 0W eee 
(9i,02, ■ ■ ■ ,&n)- The periodically continued delta 
function j(x) = J2m=-oo^( x + 27rm) accounts for 
angular periodicity, 9 = 9 + 2irm. The velocities 
yi N ) = (vi, V2, vjv), are given in terms of angle 
variables, v, = «o(cos 9i, sinflj). The collision integral 
contains integrations over the pre-collisional angles 9j 
and over N independent sources of angular noise. This 
equation is exact and can also be interpreted as the 
discrete-time analogue of the Liouville equation. Equa- 
tions of this type have been used before, for example, to 
analyze particle-based simulation methods for fluid flow 

MM- 

Assuming that particles are uncorrelated prior to col- 
lisions, the probability distribution can be expressed as 
a product of identical one-particle probability distribu- 
tions: P(9^ N \X^) = Edli-Pi Thi s approxi- 
mation of molecular chaos (MC) is valid at moderate and 
large noise strength rj and large mean free path A = r vq 
compared to the effective interaction radius R e s- It can 
be seen as a dynamic mean-field approximation because 
it neglects pre-collisional correlations. 

The assumption of large mean free path, A = 
X/R e ff 3> 1 is not very realistic for a system of swarming 
agents because it would allow agents to bypass others at 
very close distances without including them in the sub- 
sequent interaction. Nevertheless, the MC ansatz is very 
useful for several reasons. First, approximations some- 
times turn out to have a much larger range of validity 
than expected. For example, in a particle-based simula- 
tion method for fluid flow [26[, the MC approximation 
gave correct results for most transport coefficients down 
to A = 0.1. Second, the MC ansatz generates the lowest 
order terms in an expansion of the exact kinetic theory in 
the parameter e = 1 /A and thus can be seen as the first 
step towards a more complete theory. In this paper, we 
analyze these MC contributions and leave higher order 
terms for future work. 

The usual procedure 0, HU to derive a kinetic equa- 
tion for the one-particle distribution function, f(9, x, t) = 
NPi(6,x,t), is to multiply the N-particle equation with 
the microscopic one-particle density, Y^i 6(6— #j)£(x— Xj), 
for the field variables (9, x), and to integrate over all 
phases, that is all particle positions Xi and angles 9i. 
The left hand side of eq. ([5]) reduces then to NP\(9, x + 
tv,( + t) = f(9, x + r v, t + t) because integrating out 
k phases leads to N — k particle probabilities such as in 
the following example, 
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The collision term on the r.h.s. of Eq. ([5| becomes 
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FIG. 1: Illustration for the derivation of Eq. (© with M = 5 
collision partners. The selected particle 1 is fixed to the center 
of the circle of radius Rj. Particle 2 is integrated over the 
ring with inner radius Rj and outer radius Rj+i = Rj + AR 
whereas the remaining M — 2 collision partners are integrated 
over the inner circle. 



The overall prefactor N results from the fact that the 
particles are physically identical and thus every 8(9 — 
9i)5(x — term in the one-particle density yields the 
same contribution. In Eq. ([7]), without loss of generality, 
particle 1 was chosen to play a preferred role. Its position 
Xi is fixed to the field point x but the positions of the 
other particles 2, 3, ... N must be integrated over. These 
integrations are more difficult than they appear because 
the average angle $i has an implicit dependence on all 
particle positions. For example, in a system with M = 4 
interaction partners, if particles 2,5, and 7 happen to be 
the closest ones to particle 1, only they are included in the 
calculation of the average angle, $i = $1(^1,^2,^5,^7), 
but none of the others. If for example, x 2 is moved fur- 
ther away from x, another particle, say number 9, might 
become one of the closest four and will replace particle 
2 in the calculation, thus $1 = $1(^1, 9g, §5, 9j). For 
yet another spatial arrangement of particles, $j is deter- 
mined by a different set of M particles. That means, 
the 2N — 1 position and angular integrals are not inde- 
pendent of each other, they are coupled by the singular 
collision kernel 8(9 — £ — <&i). Fortunately, it is possible 
to rearrange the collision term and to reduce the number 
of integrals from infinity (in the thermodynamic limit, 
TV —> 00) to the small finite number 2M — 1. In contrast 
to other kinetic approaches plj no additional approxi- 
mations are required (see also Appendix D). 

The main idea is to draw a set of concentric circles 
with radii Rj = j AR, j = 1, 2 ... 00, around particle 1 
at xi = x. The circle distance, AR, will eventually be- 
come infinitesimal. Next, one picks M — 1 particles out 
of the available -/V particles, puts one of them in the ring 
j (between Rj+i and Rj) and distributes the other M — 2 
ones inside the circle of radius Rj, see Fig. [TJ The re- 
maining N—M particles are placed outside the ring Rj+i- 
There is (N—l)(N — 2)\/((M — 2)\(N — M)\) possibilities 
to segregate particles into three general areas: outside, 
inside the inner circle and inside the ring. The outside 
particles are allowed to move over the entire space but 
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without crossing into the ring or the inner circle. The 
single particle in the ring can move inside the ring only, 
and the inner circle particles are allowed to take any po- 
sition within that circle but cannot cross into other areas. 
For fixed ring label j these rules describe all possibilities 
to have exactly M interaction partners within a circle of 
radius Rj+i but not within a smaller circle Rj. Situa- 
tions where there is more than one particle in the ring 
are irrelevant in the limit AR — > because the probabil- 
ity for these events goes faster to zero than the one for 
single occupancy of the ring. By increasing j to j + 1, 
that is by going to the next larger ring and redistribut- 
ing particles into the three zones, one realizes that none 
of the new configurations could have already occurred at 
smaller j. We have thus constructed an alternative de- 
scription of the particle positions in terms of ring number 
j, particle labels and positions inside three distinct zones. 
For AR — > this representation is completely equivalent 
to the original description in terms of (xi,X2, . . . , Xjv), 
meaning that no double counting or omission of config- 
urations occur. This new description is crucial for the 
evaluation of the position integrals because for every con- 
figuration the identity of the M interacting particles is 
fixed. The collision integral, eq. ([7]), can now be split in 
an infinite sum over the radii Rj where M particles are 
inside Rj + AR and where N — M particles are outside, 



particles can be performed, 
1(0, *) 
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where Pj is the contribution from the outside particles, 
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Defining the average particle number in a circle of radius 
Rj+i centered around x as 



Mj(x) = / p(x')dx' 
JC(j+i) 



(14) 



and using the fact that integrating the density over the 
entire space is equal to the total particle number N, 
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p(x') dx 1 = N 



(15) 



all space 
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The index R(j) at the first spatial integral denotes inte- 
gration over a thin ring centered around x with inner and 
outer radii Rj = j AR and Rj+i = (j + 1) AR, respec- 
tively. The index C(j) means that the integration goes 
over the entire interaction circle with radius Rj , and 0(j) 
denotes integration over the area outside a circle with ra- 
dius Rj+i- The combinatorial prcfactor can be simplified 
for large N as 
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(9) 



and using the definition of the particle density p as the 
zcroth moment of /, 



f(e,x)de = p(x), 



(10) 



the integration over the angles 8m+i ■ ■ ■ @N of the outside 



one finds that 



o(j) 
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which, in the thermodynamic limit, gives 

N-M 
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Combining Eqs. (fT12|) and (fT?) leads to the final evolution 
equation for the one-particle density of the metric-free 
model 



/(0,X+Tv(0),t + T) = 



lim 
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M 



AR^O (M-2)\J_ v/2 7] J Q 

/ dx 2 / dx 3 ...dx M (18) 
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Eq. ((18)) has a highly nonlocal and nonlinear colli- 
sion term. For example, the exponent Mj is a func- 
tional of the density p which is itself a functional of /, 
Mj = Mj- [/>[/]]. However, this equation is still analyti- 
cally tractable. Note the integrations across interaction 
radii which describe collisional momentum transfer - a 
key feature of the Enskog equation - and absent in Boltz- 
mann approaches [7l [2l|. see Appendix D. 
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IV. HOMOGENEOUS SOLUTIONS AND 
PHASE DIAGRAM 

It is useful to first study homogeneous stationary solu- 
tions, /(6>,x, i) = f(9), of the Enskog-like kinetic equa- 
tion dTHJ) . The integrands are now independent of posi- 
tion, and, for infinitesimal AR, the ring integral and the 
circle integrals in (Tig)) can be replaced by 



where 



/ dx 2 ~ 2irRjAR 

Jru) 



C(j) 



dx 3 = nRj and 

Mj = n(Rj + AR) 2 p a ~ nR?p . (19) 



For AR — > 0, the sum J2j ^-R goes over to the integral 
J dR, and after substituting z = R^/npo one finds for 
the r.h.s. of Eq. ([IS]). 
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The integral over z is solvable for all M , 
2 



- 2)! J 



,2M-3 



exp(— z ) = 1 



(20) 



(21) 



and the kinetic equation (fT8|) takes the form of a nonlocal 
fixed point equation for the function f(0), 



no) = m 



M-l 

Po 



(22) 



This equation constitutes a nonlinear singular Frcdholm 
integral equation of the second kind. 



A. Disordered state 

From direct numerical simulations of VM-like models 
@, Sand from previous analytical work on the regular 
VM Q, we expect / = /o = po/(2n) to be a fixed point of 
the collision integral, independent of noise strength and 
partner number M. This constant solution describes the 
disordered phase of the model because it does not depend 
on the angle, and thus every flying direction 9 occurs 
with the same probability. In order to check this expec- 
tation we expand the collision operator of Eq. ([22]) into 
an angular Fourier series which regularizes the singular 
collision kernel 8(9 — £ — 



1(9) = Co + [CkCOs(k9) + H k sm(k9)} 



(23) 
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I(9)d9 



cos(k9)I(9)d9 



sm(k9) I(6)d6 for k > (24) 



Since /o does not depend on 9, only integrals of type 
J 21T cos(A:$i) d6\ . . . d9 M and J Q 2n sin(fc$i) d9~i . . . d9~u 
occur when / = po/(27r) is inserted into \Pm in Eq. ([20]) 
and Cfc and H k are evaluated. These integrals vanish 
for nonzero k because the average angle $i takes all val- 
ues between and 2tt with the same probability. Hence, 
Ck = Hk = for k > and it remains to find Co- After 
integrating over the noise and the pre-collisional angles 
one obtains J Q 27T d9 



Pq 1 and 



n _ Po _ f 
Co - ^ - fo 



(25) 



Thus, we indeed find that / is always a fixed point of 
the collision operator /, 



I[fo] = h 



(26) 



This is a nice confirmation that the alternative repre- 
sentation of the particle configuration in terms of rings, 
described in the previous section, is correct, and that no 
relevant configurations were left out or overcounted. 



B. Ordered state 

An ordered state of self-propelled particles is charac- 
terized by particles which have the same nonzero average 
flying direction. Such a state breaks the rotational sym- 
metry of the model and represents another fixed point 
of the integral equation ([2"2"]l . Its one- particle density, 
ford, depends on the angle and has a maximum at some 
arbitrary angle 9 which is the direction of ordered mo- 
tion. We choose 9 = because then only Fourier cosine 
coefficients are needed, 



ford(9) 

The integral equation 
algebraic equations, 
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k=0 



gk cos(k9) 
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reduces to an infinite set of 



gk = C fe (5o,3i, • ■ -ffoo) 



(28) 
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where the Ck are the Fourier coefficients of the collision 
operator 1(9), see Eq. The calculations for k = 

are identical to the ones for the disordered phase analyzed 
above, thus go = Cq = po/(2tt). To proceed, we assume 



R 



Ml 2 


3 


4 


5 


10 


Kl 1/2 1/77 


0.2624 


0.2249 


0.2008 


0.141 



TABLE I: Analytical (M = 1,2,3) and numerical results 
(M > 4) for the integrals defined in Eq. (|30[) . The asymptotic 
behavior for M ->■ oo is, ~ y/w/jteM). 



that near a specific value of the noise, r\ = rjc, only the 
lowest Fourier modes are relevant, 

.go > 9i > 92 > ■ • • (29) 

which can be easily verified a posteriori. To find this 
critical noise r]c, all terms with fc > 1 are neglected, and 
only the equation for fc = 1 in (|28|) is evaluated. Inserting 
/ = pq/(2tt) + gicos(0) into the collision operator and 
solving Eq. ([24]) for C\ yields 



.9ir + O(gi52.9 M " 2 ) 
^sin|A-(M) 

d#i . . . (i^M cos $i cos 8 1 



(30) 



1 



(27T) 



Setting Ci = gi, 77 = 77c in Eq. (|30j) and using the 
asymptotic vanishing of all higher modes, Eq. (|2TJ|) . at 
the critical point leads to an implicit equation for the 
critical noise, 



rfac) = 1 



(31) 



The Ai-dimcnsional integral K c was calculated numer- 
ically for 2 < M < 20 as well as analytically for 
M = 1,2,3 and M — > 00, see Appendix C and also Ta- 
bles U and II. These calculations are very similar to the 
ones for the network model of Ref. [l(| ■ 

The critical noise is plotted as function of partner num- 
ber M in Fig. [2] For large M, this mean field phase dia- 
gram agrees well with the one of the regular Vicsek model 
[6|, and shows the same asymptotic behavior rjc — > 2-zr 
for M — > 00. This is because for large M, the probability 
distribution for the actual radius of interaction becomes 
very narrow, thus the observed radii of interaction are 
very close to the average radius i? ff- Alternatively, this 
can be explained in the context of the regular VM. There, 
for large density, the actual particle number in a circle 
with fixed radius is very close to the average number M. 
However, even at very large M, when the critical noise is 
almost identical for both models, there remains to be an 
important difference in the stability of the ordered phase 
which we will analyze further below. As a result, the reg- 
ular and the metric-free VM become only asymptotically 
identical, at infinite M. 



C. Order parameter calculations 

The physical meaning of T, defined in Eq. (|30[) . can be 
understood by expressing the average particle momcn- 
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FIG. 2: The critical noise r\c of the metric-free VM calculated 
from Eqs. (|30l I31[) as a function of the number of collision 
partners M in comparison with results for the regular VM 
from Ref. Direct simulation results for N = 5000 particles 
and large A = \/R c h = 5.66 are shown by the full circles. If 
the noise r\ is below the symbols, the system is found in the 
ordered phase with a non-zero total momentum. The dashed 
line just serves as a guide to the eye. 



turn w as the first moment of /, 

1-2-K 

w = pu= / v(6)f(8)d8 
Jo 



v = i>o(cos#, sin#) , 



(32) 
(33) 



where u is the local macroscopic velocity. According to 
Eq. (|27|) we also have 

91 = - r C os8f(8)d8= (34) 

7T Jo T«0 

Hence, the Fourier mode g\ is proportional to the x- 
componcnt of the momentum whereas the y-componcnt 
is zero in the special case of 8 = considered here. The 
quantity T can be interpreted as the amplification fac- 
tor of momentum which can always be locally created 
or destroyed because the collision rules, Eq. (QJ, do not 
conserve momentum. One finds that T < 1 for rj > rjc- 
Thus, above the critical noise, a small nonzero momen- 
tum quickly goes to zero, and the system reaches the dis- 
ordered phase with all g^ = for fc > 0. Below the critical 
noise, rj < rjc, the amplification factor F is larger than 
unity, an initially small momentum is amplified, higher 
modes 32 , ■ • ■ <?oo get excited until a stationary state with 
nonzero momentum is reached. In order to quantitatively 
describe this ordered state and to determine whether the 
order-disorder transition is continuous or discontinuous, 
the first few members of the hierarchy of equations have 
to be analyzed for 77 < ijc ■ The smaller the noise 77, the 
more members have to be included in order to achieve 
acceptable accuracy. It is convenient to normalize the 
Fourier modes with 2go = po/ir, 



G u 



9k 
230 



(35) 
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because, for zero noise, all modes gk, k > become equal 
to 2 go. The normalized mode Gi corresponds to the order 
parameter | J2iLi v * I /( v oN) typically used in direct sim- 
ulations of flocking models @, @] ■ In order to keep track 
of the relative sizes of terms, we introduce the book keep- 
ing parameter e and assume the scaling Gk ~ e fc . This 
scaling was used previously Q and can be easily verified 
after the modes Gk have been calculated, sec Fig. [3] In- 
cluding terms up to order e 7 , the first six members of the 
fixed point equations in angular Fourier space, Eq. (]28[). 
for M — 2 interaction partners were found as, 



Gq 

Gi 
G' 2 
G, 
Gi 
G, 
A h 



1 

2 

— (2GoGi — —G\G% + 7G2G3 — -G3G4 
7r 3 5 7 

-r G i 



4 

7T 

.4 



(G1G2 — -rG()G3 4- -G1G4 — -G2G5 
3 5 7 



4 J 

— -(G2G3 — -G1G4 4- 7G0G5 — -GiGg 
7T 3 5 7 

2 M+1 /njfe 
— : — sin — 
-qk \ 2 



with 
(36) 



C. However, even the numerical evaluation is very time 
consuming for high integral dimension M > 7 and for 
high mode numbers k and i n . 

If the fixed point hierarchy, Eq. ([57)) is truncated at 
level k = kx and only terms up to order e kT are included, 
a single algebraic equation of order kx — 1 can be derived 
for the order parameter G\ . This equation was solved for 
various truncation levels, see Fig. [3J a) for M = 2, Fig. 
1151 for M = 3, and Fig. [131 for M = 7. Even for relatively 
large fcy = 5, the accuracy quickly detoriates if the noise 
is smaller than about 50 % of the critical noise, as seen 
in Fig. [3] a). In order to overcome this restriction, the 
integral equation, Eq. ([22]). was solved directly by an 
iterative numerical procedure which resolves the distri- 
bution function using 500 angular modes, details will be 
given elsewhere [27| ■ In Fig. [3] one sees that this gives 
excellent accuracy even very close to zero noise where 
G\ = 1 is predicted analytically. Once Gi is known, 
all the higher modes can be calculated and serve as the 
ground state solution in the stability analysis presented 
in the next section. 

Evaluating the algebraic equation for Gi near the crit- 
ical noise gives 



The r.h.s of Eq. (|36p contains only quadratic terms be- 
cause there is only binary interactions. For M = 3 only 
cubic terms appear because all collisions are of three- 
body type. 

The general structure of the fixed point equations for 
arbitrary M and k = 1, 2, ... 00 is given by 

00 

G k = A k J {k \hi2h---iM)G ll G l2 ...G tM 

{ij=0} 



(27T) 



M 



2k 



c / d0 {M) cos(fc$)cos(ii0i)...cos(i M 6 f M) (37) 
'0 

with <# M ) = d§! d§ 2 ... d0 M - The Kronecker symbol 
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(Af) _ 



(38) 



emphasizes that an angular integral J^(«i«2*3 ■ ■ • *m) is 
nonzero only if some addition or subtraction of its lower 
indices is equal to the hierarchy level, k = ±zi±i2 . . . im- 
For example, for M = k = 3 the term G3G2G2 would 
appear in the equation for G3 since 3 4-2 — 2 = 3 but 
not the term GfG2 because ±1 ± 1 ± 2 is never equal 
to 3 for any combination of plus and minus signs. We 
expect this property to be a consequence of rotational 
symmetry but were not able to find a mathematical proof. 
These angular integrals are M-dimcnsional and can be 
evaluated analytically for certain cases such as M = 3 but 
it is easier to determine them numerically, sec Appendix 



Gi 



vc 



1/2 



(39) 



thus, leading to the critical exponent 1/2. This is ex- 
pected for a mean field theory and was also found for the 
regular VM @. This means that at the mean- field level 
and as long as the ordered phase is stable against fluc- 
tuations near the flocking threshold, the order-disorder 
transition is continuous. One also finds the general scal- 
ing Gk ~ [ijic — r])/ric] k ^ 2 , which is confirmed in Fig. 
[3]b), Thus, the expansion parameter e can be identified 
with 



Vc - V 



(40) 



V. LINEAR STABILITY ANALYSIS 

For the regular VM Q and related metric models 
it has been shown that the homogeneous ordered state 
is always unstable to long wavelength perturbations in 
a small window rjs < t] < Vc right below the flocking 
threshold rjc- The strongest instability occurs for longi- 
tudinal perturbations where the wave vector k of a per- 
turbation is parallel to the average flying direction n of 
the homogeneous ground state. Once the angle between 
n and k increases beyond a critical angle, the system is 
stable. This linear instability explains the formation of 
large density bands in direct simulations of the VM d, [|| . 
These bands, however, have not been observed in simu- 
lations of the metric- free model [H, [l9[ . This leads to the 
obvious hypothesis that the ordered state of the metric- 
free model is stable. The question is whether it is linearly 
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FIG. 3: Order parameter Q = Gi and higher modes as a 
function of the noise for partner number M = 2. a) Ana- 
lytical calculations (dashed line) for G\ from Eq. (|36|l are 
compared with the numerical solution (solid line) of the in- 
tegral equation (|22[) and with direct simulations (circles) of 
N = 5000 particles at A = 5.66. The noise for the theoretical 
curves was rescaled by r/c,theo ~ 2.34923, whereas the noise 
for the simulation curve was rescaled by the slightly smaller 
critical noise r/c,sim = 2.2713. b) The first three modes Gi, 
G-2, and G3 are numerically determined from Eq. (I22C and 
plotted versus the distance to the threshold, e = (r/c — v)/ r lc 
in order to verify the scaling Gk ~ e k ^ 2 . The straight lines 
correspond to the exponents 1/2, 1, and 3/2, respectively . 
Since G3 is negative for 0.9 < ri/rjc < 1 the absolute value of 
G3 was plotted, which leads to the spurious dip at e » 0.1. 



or nonlinearly stable or whether it is unstable but only at 
huge wavelengths beyond the system sizes used in direct 
simulations. 

In a previous paper we had first derived the hy- 
drodynamic equations from the kinetic description by a 
Chapman-Enskog procedure and then analyzed the sta- 



bility of the hydrodynamic equations. Such a derivation 
is very tedious and involves several additional approxi- 
mations such as considering only small spatial gradients 
and assuming proximity to the flocking threshold where 
higher kinetic modes are enslaved to the lower ones. 

Here, we employ a much faster and more accurate ap- 
proach to the stability of the model. Bypassing the hy- 
drodynamic description completely, we directly impose 
spatio-temporal perturbations into the kinetic equation 
(|18[) and analyze their dynamics. Depending on the wave- 
length and the distance to the threshold, r\c ~ V: the re- 
sults of these calculations can be refined to the desired 
accuracy by increasing the number of kinetic modes to 
be included. The validity of hydrodynamic equations for 
only density and momentum is questionable anyway for 
models where momentum is not conserved because mo- 
mentum can be considered a slow variable only in special 
cases such as proximity to the threshold. Further away 
from the threshold there is no a priori justification to ne- 
glect higher kinetic modes since their relaxation rates are 
not much different from the one of the momentum. 

Introducing a small perturbation 

00 

Sf(9, x, t) = 2_. i^9n cos i n 6) + Sh„ sin (n0)] 



6g n (x,t) 
Sh n (x,t) 



Sg n e ik - x+ " t 
Sh n e 4k - x+wt , 



(41) 



of the homogeneous steady state f(6), the distribution 
function changes to / = / + 8f. The corresponding per- 
turbations of the Fourier coefficients of / are denoted as 
8g n and Sh n . For brevity we will omit the time argument 
in /, 5g n and Sh n in the following calculations. The col- 
lision operator is now spatially dependent and involves 
the following integrals 



f dx'6g n (x') = 2irRJ (kR)8g n (x) 
Jr 



dx'8g n {x') 



2nR 



J 1 (kR)8g n {x) 



(42) 



R 



which lead to 



<j> dx'f(9,x') = 2irRf + 2TTRJ (kR)5f{6,x) 



dx'/(0,x') = TTi? 2 / 



2ttR 



J^kR) <5/(0,x).(43) 



where Jo and J\ are the Bcsscl functions of the first kind. 
Note, the integrals § R dx! and J R dx.' integrate over the 
circumference of and the area inside the circle with ra- 
dius R centered around x respectively. Therefore these 
integrals are still spatially dependent. The line integra- 
tions in Eqs. (|42l |43| are related to the ring integral of 
Eq. (HI) as, 



R 



dx' = lim 

Ai?->0 



dx' 
~AR 



(44) 
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Therefore, after replacing Yl'jLa fn(j) ^ X2 *he integral 
/ °° dR § R dx' , we have 

f(6,X + TV,t + T) 

i r /2 dt 



as 



(M - 2)! J_ v/2 ri 

^...dW(«-e-$i)M,x,t) 











exp 


-TvR 2 p 









Ji(fc,i?)(5p(x) 



[27Ti?/(0 2 



A/ 

n 

z=3 



2nRJ (kR)Sf(9 2 ,x) 
, 27ri? 



-j 1 (kR)5f{e l ^) 



(45) 



Expanding the exponential function in linear order 



exp 



— irR?po — '——-—J\{kR)8p 



-ttR- pa 



k 

2ttR 



J 1 (kR)Sp) +0(5p 2 ) (46) 



and integrating over the collision radius R, we arrive at 



/(6>,X + TV,t + T 

1 i rf/a 



M-l „ 



d0! . . . d6 M S(9 - £ - $0^(0!, 02, • ■ ■ hi) (47) 



where 



J~(9i,02, ■ ■ ■ 


9 m) 






= A fl /3 • 


■ Jm 


1 - (M - 


-1)^ 
5o 


+<*/i /a h ■ 


■ ■ fid 






+/l 5/2 /S ■ 


■ ■ fa 


iPi(M- 


1,1, z) 


+/i /a <5/s ■ 


■ ■ /m 


iPi(M- 


l,2,z) 


+ ... 

+h hh-- 


■Sf M 


iPi(M- 


l,2,z) 


+0(S 2 ). 









(48) 



The abbreviations fj and Sfj stand for f(9j) and 
Sf(9j , x) respectively, and 1F1 (a, 6, z) is the confluent hy- 
pergeometric function with argument z = — fc 2 /(47rpo), 
see Appendix A. The nonlinear perturbations (denoted 
as 0(5 2 )) will be neglected in the following. Since the 
collision integral (|47|) is symmetric under permutation of 
the pre-collision angles 9i, the integrand can be written 




l-(M-l)-^ 1 F 1 (M, 2, z) 
5o 



l + iFi(Af-l,l,z) 



+(M-2) !Fi(A/-l,2,z) 



(49) 



Note, the integrand Eq. (|4"9")l is the general form for M- 
particles collision. For the special case, M = 2, where 
there is no particle inside the inner circle, this simplifies 

5 go 



to F(9 U 9 2 ) = fi f 2 [1-^ e z J + 5h / 2 (1 + e z ), where 
iFi(a, a, z) = e z has been applied. 

A. Fourier expansion of the collision integral 

The strategy to derive the growth rate ui(k) is to ex- 
press both sides of Eq. (|47j) in terms of the angular 
Fourier coefficients C„(x, t) and H n (x,t), in analogy to 
the homogeneous case, Eq. (f2~4")) . By equating the dif- 
ferent expressions for C n and H n obtained from the left 
hand and the right hand side of Eq. (|47[) . a matrix equa- 
tion for the perturbations Sg k and Sh k is constructed. 
The dispersion relation w(k) follows from demanding 
that there is a nontrivial solution. 

To check the validity of our approach, in particular 
expansion Eq. (|4"6"|) . we first calculate the zero mode Co 
of the collision integral on the r.h.s. of Eq. (|47|) . namely 



Cn 



1 

2^ 



i[f(9,*,t)} de 



(50) 



= <?o + <55o[l + iFi(M-l,l,z) 

+ (M - 2) !Fi(Af - 1, 2, z) - (M - 1) iFj (M, 2, z)] 

The confluent hypergeometric functions cancel according 
to the identities 

1 F 1 (a,b-l,z)- 1 F 1 (a + l,b,z) = (51) 

x iFi(a+ 1,6+ l,z) 

iFi(a + 1,6, z) - iFi(a,&,z) = ^ iFi(a + 1, 6 + 1, z), 

b 

and one finds 

C (x,t) = g + 5g (x,t) = ^p(x,t) (52) 

Z7T 

This is the expected result because the collisions do 
change the local velocity and other moments but they 
do not modify the local density p. The local density is 
the zeroth moment of /. Thus, the invariance of density 
under collisions requires 



p = J f(9)d9 = J I[f(9)]d9 = 2irC 



(53) 
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which is exactly what we found in Eq. (|52[) . 

The non-zero modes of the collision integral are 

C n ^o = A n J d6i...de M J r (0i,02,---§M)cos(n^), 

H n ^ = A n [ dh---de M T{hA,---(>M)s™{n§). 



(54) 



where 



^» = -Ti^T— ~ sin (-pl ■ ( 55 ) 
p™ 1 1 nnri V 2 / 

For a homogeneous ordered state close to the flocking 
threshold, one has g n ~ e™ ~ I 2^- T> J ; where e mea- 
sures the relative distance to the critical point, see Eq. 
(|40jl . Suppose we want to expand C n and H n up to order 
of e z , we will need all the integrals 



ddx . . . d9 M T {n<S>)T 1 {k 1 9 1 ) . . . T M {k M e M ) (56) 



which satisfy ki + fc 2 + • • • + fejvf < 2, where fcj is a 
non-negative integer and function T,(x) is either sin(cc) 
or cos(x). However, not all the integrals have non-zero 
value. We found that integrals which do not satisfy the 
condition, ±n ± ki ± fe 2 i • • • ± = 0, vanish. Fur- 
thermore, if the total number of sin function inside the 
integral, including sin(?i<I>), is odd, the integral also van- 
ishes. For the binary collision case, M = 2, by defining 

K%£ q = (cos (n$) cos (p9i) cos {q6 2 )) 

K™ s q ee (cos (n$) sin (pflO sin (<?0 2 )) (57) 

K^ s g ee (sin cos (p^) sin (g0 2 )) 

with (. . .) ee dOi J Q 27r c?6*2 / (27r) 2 , we have for odd n, 



TfCCC 

npq 

TSCSS 

npq 

TSSCS 

npq 



a 1 



(+01 +a 2 + a 3 )/(27r) 
(-01 + a 2 + a 3 )/(27r) 



$p-\-q,n 



"2 



0-3 



(+ai 






-1 


P - 







-1 


P + 


— 




-1 


P + 


— 

'7 



$p—q,n 



$q~p,n j 



(58) 



where i is the imaginary unit. When n is even we find 



If ess 
npq 



ry-scs 
npq 



-^T — ^2p,7i^2q,7i 



(59) 



B. Fourier expansion on the left-hand side 

So far we have considered the Fourier expansion of 
the collision integral, which is the right-hand side of the 
Enskog-likc equation (|4"T|) . On the left-hand side, writing 
down the Taylor expansion around (x,£), we have 



f(9, X + TV, t + T) = V — (d t + V a d a ) n f{9, X, t) 
n=0 

00 

/ + e T ^+^ v) £ \8g q cos (q9) + 6h q sin (60) 



9=0 



The wave vector k is split into a longitudinal part and a 
transversal part with respect to the average direction of 
the ordered state n, k = fciin+fcxt. Since the coordinate 
system was chosen with normal direction n = (1,0) and 
transversal direction t = (0, 1), one finds 



k • v = vq (fci I cos 9 + fcj_ sin ( 



(61) 



The identities 



J (x) + 2 ]T J n (x)i n cos (n6) (62) 

n=l 
00 

Jo(x) + 2i 22 J2n-i{x) sin (nO) 

n=l 

00 

2 ^ J 2m {x) cos (m9) , 



(63) 



are used to express the factor e rlk v in Eq. (|60|) in terms 
of Bessel functions. Because of simplicity and because 
the strongest instabilities of the original Vicsek model 
occur in the longitudinal direction, where k = feiin, we 
restrict ourselves to this case. It is straightforward to 
generalize the following analysis to arbitrary directions 
of the wave vector. Using Eq. (f6"2"| we rewrite Eq. ([S0|) 
as 



f(9,x + rv 7 t + r) = ^2 gj cos jt 



3=0 



Joirvk) + 2 V] J p (Tvk)i p cos (p0) 



P =i 



00 

^ cos (q#) + sin (gfl) 



9=0 



(64) 



This can be further converted into the Fourier series 

/(0,x + TV,i + r) = C + Cicos6i + C2cos(26i) + ... 

+ Hi sin 6» + H 2 sin (2(9) + . . . , (65) 



11 



and allows us to read off the coefficients C n and H n . 
Since, for non-negative integers p, q, and n 



(cos (p9) cos (q6) cos (n9)) 



(cos (p0) sin (qO) sin (n9)) 



Sp-\-q 7 n &p—q,n ~t~ &q—p,n 



Sp-\-q 7 n &p—q,n ~t~ ^q—Py 



(66) 



with (...) = J Q 27r d6/ '(2tt) we have 



Co 



1 f 2 " 

— y f(9, X + Tv,t + r) d9 

CO 

<7o + e™X>%<^ 



p=0 



1 f 2 " 

- / f(0, x + tv, i + r) cos (n0) d6> 
Jo 

fl n + e™ [ Jo Sg n (67) 



^ ^ ft 'Jp $Qq ($p-\-q,n &p—q,n H~ ^— p+<?,n) 

CXD 

9n + e™ £ { iln ' qlj \n- q \ + * n+q Jn +q ) Sg q 

q=Q 

1 

-/ f(6, x + rv, t + r) sin (n<9) d(9 
Jo 

e™ [ Jo <5/i„ 

OO 

^ ^ ft Jp Shq {&p-+-q.n ~ ^p—q.n ~\~ & —pA-q^n) 
p=l,g=0 



= e™ ^ J|„_ g| - ? "+«J„ +f; ) fli,. (68) 

In the following, we show examples of C„ and H n in the 
expansion up to order e 2 , 



Co = 
Ci = 

+ 

c 2 = 
+ 

Hi = 
H 2 = 



go + e TUJ [J Sg + iJiSgi - J 2 Sg 2 + 

gi + e™[2iJiSg + (J - J 2 )5gi 

i(Ji - Js)Sg 2 + ■■■} 

g 2 + e™[-2J 2 5go + i(Ji-J 3 )5 gi 

(Jo + Ji)dg 2 + ...} 

e™[(J + J 2 )5h l + i{J x + J 3 )5h 2 + 

e TU [i(J! + J 3 )Sh 1 + (J - J 4 )5h 2 + 



(69) 



Equating the expansion of the left-hand side, Eqs. ([6"TI 
, to the one of the right-hand side, Eq. (|54[) . a matrix 



equation of the following general structure is found, 



( Coo Coi • • • Aoi Aq 2 
Cio Cn • • • An Ai 2 



Bio Bn 
B 2 o Bn 



H\\ H X2 
H 2 i H 22 



\ 



Shi 
Sh 2 

J V : J 



(70) 



According to Eq. for the case k = k|| = kx con- 

sidered here, the Sgi represent longitudinal perturbations 
which describe a change of the magnitude but typically 
not of the direction of mean flow. The Shi stand for 
transversal perturbations which describe odd variations 
of the distribution function, 5f(6) = —Sf(—6), and mod- 
ify the flow direction. Since we only consider linear per- 
turbations, the C n of Eq. (|67j) are found to only depend 
on Sgj but not on Shj. Similarly, H n depends only on 
Shj. Therefore, the block matrices A a p and B a p are zero. 
That means, the Sgj, are decoupled from the Shj. 

The modes 5g n and Sh n are neglected for n > nc and 
the block matrices C a p and H a p are truncated corre- 
spondingly. This truncation is motivated by the obser- 
vation that, for any nonzero noise, the angular Fourier 
modes, <?„, of the homogeneous ordered state decay to 
zero with increasing mode number n. It is plausible to 
assume that the perturbations of these modes, Sg n and 
Sh n , show a similar behavior. Of course, this decay will 
be quite slow if the noise is much smaller than the criti- 
cal noise r\c ■ This requires a sufficiently large truncation 
mode number nc- Its correct choice is discussed in the 
following chapter and shown in Fig. ITT1 

Setting the determinants of both matrices equal to zero 
leads to 2nc — 1 different branches of the dispersion re- 
lation w(k). The real part of the growth rate, Rc(oj), of 
a few of these branches is plotted in Figs. [4] -Qj] for dif- 
ferent distances to the threshold and for various partner 
numbers M . 



C. Results and discussion 

In order to analyze the dispersion relation we distin- 
guish between longitudinal and transversal modes. The 
longitudinal modes are shown in the left panels of Figs. 
S] - [8J and [10] and are characterized by changes in the 
density Sp tx Sgo and in the ^-component of the average 
flow Sw x oc Sgi. These changes are always accompanied 
by corresponding perturbations of the higher order an- 
gular coefficients, Sg 2 , Sg 3l . . . Sg nc -\. The right panels 
depict the growth rates of the transversal modes. We 
further distinguish between hydrodynamic and kinetic 
modes. The hydrodynamic modes correspond to exactly 
conserved quantities and go to zero for k — > 0. In Figs. 
141-151 and [TU1 thev are plotted as solid lines. The kinetic 
modes, defined as those going to a non-zero value at zero 
wave number, are given as dashed lines. The figures show 
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that there are only two hydrodynamic modes. This is ex- 
pected because only two quantities are exactly conserved 
in the collisions: mass and the individual kinetic energy 
of every particle since the particles' speed never changes. 
A closer look reveals that the hydrodynamic mode de- 
picted in the left panels of the figures is a sound mode 
which is related to longitudinal changes of x-momentum 
and density. The hydrodynamic mode in the right panels 
of Figs. l4l-[T0lcorrcsponds to changes of the y-component 
of the momentum. At zero k this mode describes the re- 
sponse to a small rotation of all particle velocity vectors 
by the same amount. Due to rotational invariance this 
is a Goldstone-mode which meets no resistance and lu is 
zero. 

We label the top dashed line in the left panels as a 
"pseudo" -hydrodynamic mode because to(k = 0) is zero 
only at the critical point, r\ = r\c but is negative away 
from this point. This mode is related to the fact that, in 
general, momentum is not conserved in VM-like models 
but at the critical point there is no amplification of mo- 
mentum perturbations. Another way to understand this 
is to look at the hydrodynamic equation for the momen- 
tum of the VM, Eq. (5) in Ref. Q, which has the general 
shape 

d t w = (r - l)w + 0(w 2 w) + . . . (71) 

where Y is the amplification factor defined in Eq. (|30p. 
At the critical point, Y(rjc) = 1 and w is small. Thus, 
at rj = rjc, in linear order in the momentum density w, 
momentum is conserved. 

In the left panels of Figs. 0] - [H we observe that the 
"pseudo" -hydrodynamic curve drops to lower negative 
values the further away one is from the critical point. 
In Fig. [6j at r\ = 0.6?yc, we are so far away from the 
critical point that this mode has now similar relaxation 
rates as the other purely kinetic modes. 

The main result of this linear stability analysis is that 
there is no longwave instability. In particular, we found 
that Re(u>) of the sound mode, like all other modes, is al- 
ways negative at small wave numbers. This is in contrast 
to the regular VM whose modes we show for comparison 
in Fig. [TU] For the VM, the sound mode is clearly unsta- 
ble at wave numbers below fee, confirming the previous 
result from Ref. @ which was based on hydrodynamic 
equations. In the right panels of Fig. [lOlwe see that all 
transversal modes are stable. More details on the regular 
VM will be reported elsewhere, [27| . 

For small truncation level nc we did see an instability 
at higher wavenumbers, for kX gl 2.6. However, when nc 
is increased, this region is shifted to even higher k 3> 1 
whereas ui remained unchanged at low wavenumbers, see 
Fig. HU This strongly suggests that the short wavelength 
instability is spurious. It is just a result of neglecting 
higher order terms in Eqs. (|37[) and/or |70p that can 
be easily remedified. This result is also consistent with 
direct simulations of Ref. fl9j which showed no sign of 
instabilities at any k. For partner number M = 2 we in- 
vestigated how far from the threshold one can go and still 



have linear stability of all modes. Figs. [U [5] and [6] show 
that no instability occurs down to n = 0.6r]c- We did 
not go to even lower noise because many more terms of 
the ground state solution G n and also more perturbation 
modes 8g n would be needed to achieve reliable results. 
Finally, we were interested in how the large partner num- 
bers seen in experiments 0, [22j modify linear stability. 
Calculating w(k) for all M between two and seven, see 
Figs. [3 [8] and [H it is clear that the situation becomes 
even better: the larger M is, the more stable the modes, 
especially the sound mode, become. 

To conclude, at M = 2 - 7 there is no signs of linear in- 
stability in the metric-free model at or below the flocking 
threshold. This supports previous claims based on direct 
simulations [lj| that the order-disorder transition in this 
model is continuous and is not made discontinuous by 
linear instabilities near the threshold. 

To understand why the long wave-length instability of 
the regular VM does not appear in the metric-free ver- 
sion, let us compare the momentum amplification factor 
T for both models. According to Eq. (f30j). for the metric- 
free case, r depends on the partner number M which is 
a constant. Thus T is the same in regions of low and 
high particle density. This is not the case for the met- 
ric VM. The amplification factor (which is defined as A 
in Eq. (3) of Ref. 0) depends on the local number of 
collision partners Mjj(x) which is proportional to the lo- 
cal density. Analyzing this expression shows that, in the 
metric case, Y is monotonically increasing with density 
and at high density scales as Y ~ J~p. As a result the 
critical noise rjc also increases with density. A possible 
explanation for the longitudinal instability of the regu- 
lar VM goes then as follows: Assume a spatial region of 
low density, p^ < p a . The local critical noise, tjc(pl) 
is small, and hence this region corresponds to a point in 
cither the disordered part of the phase diagram or to a 
point which is only slightly below the flocking threshold. 
This means, on average, particles go in almost all direc- 
tions and the macroscopic local velocity w/p is small or 
zero. This is consistent with a small amplification rate 
r«l and also with Figs. [3^) and [131 which show that 
the average speed decreases monotonically if the flocking 
threshold is approached from inside the ordered phase. 
In a dense region with pu > pa, the local critical noise, 
t]c{pd) would be significantly larger than n. This region 
is then described by a point deep inside the ordered part 
of the phase diagram. Thus, particles would be strongly 
aligned and the average local speed would be large, con- 
sistent with a large Y. These are exactly the conditions to 
form a density wave: particles in high density regions are 
more aligned and "invade" regions of lower density where 
particles perform a slower average motion and thus are 
"not organized enough" to escape from the dense crowd 
coming in. In the metric-free model, the local density 
does not couple to the average particle motion in this 
way. The mechanism for density wave formation dis- 
cussed above is absent. A similar discussion of the ab- 
sence of instabilities in metric-free models is given in Ref. 
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FIG. 4: Real part of the growth rate ui as a function of 
wave number for M = 2, close to the flocking threshold, at 
■q = 0.997/c, calculated from Eq. (|70|l . Other parameters: 
A = 2, po — 0.1. The ground state solution is accurate up to 
order e 5 . Perturbations Sg„, Sh„ with n > 6 were neglected. 
Part a) shows solution branches for the determinant equa- 
tion, det(C) = 0, where the block matrix C a p is defined in 
Eq. (|70[) . These curves describe different longitudinal exci- 
tations where typically all angular perturbation coefficients 
8go,Sgi, . . . Sg n are nonzero. Within a particular excitation 
mode and for a given wave number, the coefficients Sg n occur 
in fixed specific ratios to each other. Part b) shows the real 
part of the growth rate for different transversal modes which 
are composed of the 5hi,5h2, . . .5h n . Hydrodynamic modes 
are plotted as solid lines, the kinetic modes are dashed. Parts 
c) and d) are just zoomed in versions of a) and b), respectively, 
to better show the small k behavior. 



21[ . In appendix D wc relate our approach to the work 
of Ref. [2l| and investigate the role of collisional momen- 
tum transfer by considering various limits of the general 
Enskog-like kinetic equation, Eq. (JT5J). 



VI. DIRECT NUMERICAL SIMULATION 

In order to verify our analytical results we also per- 
formed direct numerical simulations of the model defined 
in Eq. ([1]). A quadratic simulation box of size Lx L with 
periodic boundary conditions is used and seeded with JV 
particles. Their initial positions and flying directions are 
chosen at random. For every particle i, the distances to 
all other particles are measured and the M — 1 particles 
with the smallest distances are defined as the neighbors 
of particle i. This differs from the neighboorhood defini- 
tion of Ref. [lj| by means of a Voronoi-construction. 



FIG. 5: Real part of the growth rate ui as a function of wave 
number for M = 2 and r\ = 0.9077c- Notation and other 
parameters are the same as in Fig. [4] 
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FIG. 6: Real part of the growth rate oj as a function of wave 
number for M = 2 far from the flocking threshold at r\ — 
0.6077c- Notation and other parameters are the same as in 
Fig. H 




FIG. 7: Real part of the growth rate oj as a function of wave 
number for M = 3 close to the flocking threshold, at 77 = 
0.9977c- Other parameters are the same as in Fig. [4] 
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FIG. 8: Real part of the growth rate ui as a function of wave 
number for M = 7 at 17 = O.QQr/c- Other parameters are the 
same as in Fig. [4] 



After the system had relaxed into a stationary state, 
we performed a time average of the order parameter 



1 



Nvn 



N 



(72) 



and plotted it as a function of various parameters, see 
Figs. [T3landfT4l While typical particle numbers are N = 
1000 and 5000, we first performed a test at N = M = 
2 because an exact result without the molecular chaos 
approximation can be derived for the order parameter 




/v 2 



if 77 < 7T 

if 7r < rj < 2tt 



(73) 



in Appendix B. As seen in Fig. [T2lthe 
N = 2 are in excellent agreement with 



see Eqs. 
simulations for 
this formula. 

Another analytical result was obtained for maximum 
noise strength 77 = 27r where a mapping to a random walk 
can be utilized and for large N 



2^) - 7 -^= 



(74) 



is obtained, see Appendix B. Our simulations were also 
in excellent agreement with this expression. The deriva- 
tion of the phase diagram relied on the approximation of 
molecular chaos. For nonzero noise and finite M , this ap- 
proximation becomes exact at infinite mean free path A. 
Therefore, the order parameter and the phase diagram 
were measured at a large ratio A = A/i?, c ff = 5.6 and 
compared with analytical results in Figs. [T3] and 1141 In 
order to investigate the importance of the mean free path 
(mfp), the critical noise at fixed M but different A was 
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FIG. 10: Regular Vicsek model: Real part of the growth rate 
u) as a function of wave number for low density, (M) = 0.1, 
close to the flocking threshold, at n = 0.997/c- Note the long 
wavelength instability in part c). Other parameters: A = 
A/ R = 4. Details will be given elsewhere [27l |. 




FIG. 11: Effects of truncating the matrix equation, Eq. (|70|l : 
Real part of the growth rate tj of the slowest Sg n mode as a 
function of wave number for truncation levels nc = 3 (dashed 
line), nc = 4 (dotted) and nc = 5 (solid line). Parameters: 
M — 2, n — 0.9977c, A = 2, p = 0.1. 



determined, see Fig. [T5l We found that for both M = 2 
and M = 7 the influence of the mfp is only negligible if 
A is above one. For smaller mfp's, the differences are sig- 
nificant and are subject to current studies. For example, 
r\c drops by almost a factor of three if A is reduced from 
5.6 to 0.1. 
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FIG. 12: The order parameter G\ as a function of noise for 
the special case N — M = 2. The exact solution, Eq. (|73|) . 
(solid line) perfectly agrees with the simulation (circles). 
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FIG. 13: Direct simulations: The order parameter Gi as a 
function of noise for various M = 2 ... 7 and at N = 5000, 
A = X/R c a — 5.66. The analytical solution for M — 3 (filled 
circles) was obtained from Eq. Q37p with terms up to order 



The main assumption of our mean-field theory, the 
molecular chaos approximation, predicts that the particle 
number in a given box is Poisson-distributed @. Hence, 
the probability to find n particles in a box of area V is 
given by 



Pn 



-<»>. 



(75) 



where (n) = Vpo is the average particle number in that 
box. In order to indirectly test the validity of the molec- 
ular chaos approximation, we divided our simulation do- 
main into 25 x 25 quadratic boxes such that (n) = 8, and 
recorded how often a box was occupied by a given particle 
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FIG. 14: Direct simulations: The order parameter for dif- 
ferent particle numbers: N = 1000 and 5000; A = 5.66. The 
theoretical result for M — 7 calculated at order e from Eq. 
(|37l) is shown by the filled circles. 
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FIG. 15: Order parameter for M — 7 at different mean free 
paths, A = X/R c ff, with TV = 1000. The lines just serve as 
guides to the eye. 



number n. These measurements were averaged over all 
boxes and over time. The results for M = 2 are shown in 
Fig. [16] and compared to an analytic continuation of Eq. 
(|75|) to real numbers where n\ was replaced by the gamma 
function, T(n + 1). For large noise rj = 6, which is far 
beyond the order-transition threshold r\c,sim = 2.27. the 
histogram shows perfect agreement with the Poisson dis- 
tribution. This confirms our expectation that the Molec- 
ular Chaos approximation should be valid in the disor- 
dered phase. However, at noise rj = 2.2, that is in the 
ordered phase but only 3% below the threshold, a clear 
deviation from the Poisson distribution occurs and the 
maximum of the curve lies about 20% below the Poisson 
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FIG. 16: Histogram for different noises at M = 2 in com- 
parison with the Poisson distribution (solid line) as a func- 
tion of particle number per box, n. Simulation parameters: 
TV = 5000, A = 5.66, (n) = 8, time average over 25000 mea- 
surements. The dashed and dotted lines are only guides to 
the eye. 



curve. This is interesting because despite this deviation 
the shape of the order parameter curve and the critical 
noise value differ much less from the mean field predic- 
tions. Finally, at very low noise, r) = 0.3 = 0.13r)c, S im, 
the particle number distribution is much wider than the 
Poisson distribution. In particular, it is much more prob- 
able to find empty boxes and boxes occupied with more 
than three times the average number. We see that even 
though no density bands are observed as in the regular 
VM, there is still anomalously large density fluctuations. 



VII. CONCLUSION 

We have presented a detailed, systematic derivation of 
a kinetic theory for a model of self-propelled particles 
with metric-free interactions. This discrete-time model 
has genuine multi-body interactions and was introduced 
in Ref. Q. The sole approximation in the derivation 
was the assumption of Molecular chaos which we used to 
reduce an exact Master-equation for all particles to an 
equation for the one-particle density. 

This novel Enskog-type kinetic equation, Eq. (jlgf . is 
one of the main results of this paper. Using this equation, 
the transition from a disordered state to a homogeneous 
state of collective motion was studied for various num- 
bers of interaction partners, M. We calculated the phase 
diagram and the order parameter analytically as well as 
numerically and also performed direct simulations. We 
found excellent agreement within a few percent between 
theory and simulation as long as the mean free path is at 
least several times larger than the effective interaction ra- 
dius. In order to test the validity of the molecular chaos 
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approximation, we recorded particle density histograms 
and measured how lowering the mean free path affects 
the phase diagram. 

Simulations of the regular Vicsek-model showed that 
the flocking transition becomes discontinuous once the 
system size is beyond a certain critical length [f|. This 
observation has been linked to a long wave length in- 
stability of the ordered phase right below the nocking 
threshold @, 0|- In this paper, we have performed a lin- 
ear stability analysis of the ordered state of the metric- 
free model in order to investigate the nature of the flock- 
ing transition in the presence of topological interactions. 
This was done by directly imposing perturbations into 
the kinetic equation without first deriving hydrodynamic 
equations. Such derivations are very tedious for models 
with a finite time step and multi-body interactions, see 
Refs. 0, H3, If one is only interested in the lin- 

ear stability of a certain phase it is more convinient to 
use the kinetic equations directly. An additional advan- 
tage is that by not imposing closure of the kinetic equa- 
tions at a predetermined low level and refraining from 
gradient expansions of any kind, higher accuracy and a 
larger range of validity of the stability analysis can be 
achieved. The main result is that for all partner num- 
bers 2 < M < 7 we tested, all modes arc stable right 
next to the nocking threshold. For select M we verified 
that even very far from threshold, no linear instabilities 
occur. This result is consistent with direct simulations 
of the metric-free model where no high density bands - 
a sign of instability - were observed and where the flock- 
ing transition was found to be continuous. While our 
results do not come as a big surprise, they do rule out 
the possibility of a linearly unstable but nonlincarly sta- 
ble ordered state. This would lead to an inhomogeneous 
ordered state which would be hard to identify in a direct 
simulation if the inhomogeneity is small. 

The existence of a longitudinal instability can be re- 
lated to the different parameters that span the phase di- 
agram: In models with metric interactions, the critical 
noise depends on local density. Thus, different spatial 
regions can be characterized by different phase space dis- 
tances, Ar] = 77c(pi oca i) — 7]. to the flocking threshold. In 
metric-free models all regions are at the same distance to 
the threshold. With the additional facts that the order 
parameter is monotonically increasing with A77 and that 
the critical noise of the metric VM is increasing with 
density, an intuitive understanding of the occurence of 
density waves in the metric models and of their absence 
in topological models can be obtained. 

It remains an open question how to systematically go 
beyond the approximation of molecular chaos in order 
to improve the results at low mfp. One possibility is 
to derive additional noise terms in the hydrodynamic or 
kinetic equations along the lines of Refs. |28l - l3l| . Work 
in this direction is in progress. 

We did not derive hydrodynamic equations for the 
metric-free model because it is fairly obvious that the 
must have exactly the same shape as Eq. (5) of Ref. 



which was derived for the regular Vicsek model. This 
is because both models have the same symmetries (ro- 
tational and translational) , no Galilean invariance, and 
the same set of conserved quantities (mass and kinetic 
energy). The Chapman-Enskog expansion of Ref. [f|, 
which is basically a gradient expansion, gave all possible 
terms allowed by the symmetries and by the order of the 
expansion. Therefore, the metric-free model has no other 
"choice" than picking the same terms just with different 
coefficients. 

We also identified the ad hoc Boltzmann-like collision 
integral of Ref. plj as the zero wave number or infinite A 
limit of our theory for the special case of M = 2, see Ap- 
pendix D. This was used to quantify the relevance of "col- 
lisional momentum transfer" to the linear stability of the 
ordered phase. Finally, the kinetic formalism presented 
in this paper might also be useful to treat other "exotic" , 
multi-body interaction rules which are often postulated 
in ecological modelling of animals [32| . human crowds 
[33l ] , and interacting robots [34l-l36j . 
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Appendix A: Hypergeometric functions 



Using the expansion of Bessel function of the first kind 



^ «!(„ + «)! V2 



and the integral 



x m e -ax 2 dx = r[(m + l)/2] ^ a > Q ^ m > 



2 a ( m+1 )/2 



we have 



(77) 



x m J n (x)e- ax dx 



(-I) 8 _L / m + n + 2s + l \ _™ ±2 ±2 £± i 

~ f^o s!(n + s)!2"+ 2s+1 V 2 

Setting I = (m + n + l)/2, the series becomes 



1 (l + S-l)(l + S-2).. .1(1-1)1 (-1 
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(78) 
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where 



1 F 1 (a,b,z) = J2 



(a) a z 8 

(b) s s\ 



(79) 



and the integration area has to be split into several do- 
mains in order to correctly treat the absolute value in 
Eq. flS3J|. This gives 



is the confluent hypergeometric function in the notation 
of Pochhammcr symbol (x) n — x(x + 1) . . . (x + n — 1). 
If the argument of Bessel function rescales as x — > kx, it 
is easy to show that the integral becomes 



x m J n {kx)e- ax dx 



1 



p ^ m-\-n-\-l ' 



2 n+l a (ro+n+l)/2 r(n + l) 
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Appendix B: Exact solutions 

A. Case N=2 

For M = N every particle is neighbor to every other 
particle, which can be analytically exploited. We con- 
sider the special case N = M = 2. The order parameter 
is expressed as, 



O 



1 



(|ni + n 2 | 



(81) 



where are the normalized velocity vectors after a col- 
lision. The angular brackets denote the average over the 
two uncorrelated angular noises £i and £2- The average 
angle $ is the same for both particles and is therefore 
irrelevant for the order parameter. We choose it to be 
zero. In this case, fij = (cos sin and we have 



1 



-77/2 J-n/2 



J(£i>60 = -^(cosa + cos£ 2 ) 2 + (sin6 + sin&) 2 



-^2(l + cos(6-6) 

„ 6 - £3 



6 6 

cos — cos — 
2 2 



.6.6 
sin — sin — 
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(82) 



For r\ < tt the integrations are straightforward since |£i 
^2! < 7T and give 
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For it < r) < 2tt new variables are introduced, 
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(86) 



with 



Ai = 
A 2 = 
A3 = 
A 4 = 



n/2 
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where c = cos(ai/2). 

Integrating over ct2 yields, 



n = - 
'I 



-■n/2 



(—c)(r] + ai)dai+ / (+c)(>? + ai) dai 
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After integration we obtain for it < r\ < 2ir: 



(87) 



1 + 7/ — 7T + COS ■ 



At 77 = 7r both expressions, Eq. (|83|) and (|88|) match at 
f2 = 8/7r 2 w 0.81. This point, rj = tt, is also the turning 
point of the order parameter curve, Fig. 1121 

At the largest noise value, 77 = 2tt, we find fl = 2 /it w 
0.6366. It is interesting to note that the approximation 
for large N and rj = 2ir, Eq. (|M|) . even works well for this 
N = 2 case, since 7/(8^2) = 0.6187 is only 3% smaller 
than the exact result. 
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B. Case N > 1, 7? = 2tt 



The order parameter formula, Eq. (|8Tj) is generalized 
to N particles, 



n = 



N 




(89) 



Using cos 2 & + sin 2 £j = 1 the terms inside the square 
root can be reordered with the result, 




(90) 



Since c,; = cos^ and Si = sin^ vary between —1 and 
1 the terms A and B will be smaller than N for most 
realizations when the particle number N is large. We 
therefore attempt a Taylor expansion of the square root 
in Eq. (|90|) whose validity can be checked a posteriori. 
We obtain. 
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(91) 



Since the angular noises £j are uncorrelated one finds 
{A) = (B) = {AB) = 

(A 2 ) = ^) . _<^_1> (92, 

Substituting into Eq. (|9"Tj) gives 

1 



1 / 1 
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(93) 



Thus, we arrive at the following approximative expression 
for the order parameter at N 3> 1: 



7 1 

fi(rj = 2tt) rj = 



(94) 



Appendix C: Integrals for M = 2 and M = 3 

In order to calculate the angular integral Kq(M), see 
Eq. (|30[) . and integrals of similar type, the average angle 



<f> is expressed by means of the local order parameter 
vector Ljv/ = (Lm x , Lm.y), defined as, 



M 

£ 

i=i 



(95) 



where fij = (cos a,, sin a,) = v^/vo is the normalized 
velocity vector for agent i. The sine and the cosine of 
the average angle follow as cos$ = L X /\L\ and sin<E> = 
L y /\L\. The average angle is given by $ = &t&n(L y /L x ). 
For M = 2, trigonometric addition rules can be used to 
simplify the integrations, 

ai + a 2 ai - a 2 
L x = cos a\ + cos a 2 = 2 cos cos 

. Ot± + 0>2 — a 2 

L. y = sm ai + sin a 2 = 2 sm cos (96) 



yielding 



Q'l + Q'2 

2 

Q'l + Q'2 



for \ai — a 2 | < 7r 

h 7T for |CK1 — CK2I > 7T (97) 



for < ati < 27r. The integral over ai and a 2 is split into 
four parts, 



27T p27T 
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Jo 

dai I / da-i 
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da 2 
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(98) 



where in the first and third part \cti — a 2 | < 7T, and in 
the second and fourth term one has \a± — a 2 | > All 
functions under the integral are now products of sine and 
cosine with a linear combination aai +ba 2 - Therefore, all 
integrals of the kind shown in Eq. (f3"0")) can be avaluated 
analytically for M = 2. For example, one finds Kq(M = 
2) = 1 /%, see Table I. More details and information about 
how to exactly evaluate collision integrals for M = 3 and 
M —> 00 will be given elsewhere [27} ■ 

In order to perform the order parameter calculations 
and the linear stability analysis for systems with three- 
body interactions, M = 3, the following integrals are 
needed: 



T^CCCC 
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TfSCSC 
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(cos(7Ti<f>) cos(p6*i) cos(g6* 2 ) cos(r6>3)) 
(sin(m$) cos(p£?i) cos(<7#2) sin(r#3)) 
(sin(m$) cos(p#i) sin(g# 2 ) cos(r#3)) 



K^ r = (sin(m$)sin(^ 1 )cos(^2)cos(r0 3 )> (99) 
with (...) = j 271 d§if^ d8 2 f^ d8 3 /(2ir) 3 . Using the 
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definitions 



Appendix D: Effect of collisional momentum transfer 



TSCCCC 


= SiK 


rnpqr 


jy-sccs 
rnpqr 


= S 2 K 


ry-scsc 
rnpqr 


= S 3 K 


ry-sscc 
rnpqr 


EE S 4 K, 



(100) 



where K and Si depend on the quadrupel (m,p,q,r), 
the numerical values of these integrals for in < 5 and 
p + q + r < 6 can be constructed from Table II. 



m p q r 



K 



Si S2 S3 S4 









1 

2 

1 1 
1 1 
1 2 

1 
1 

2 

1 1 
1 1 

1 2 

2 2 

1 

1 

1 1 
1 1 
1 2 

1 





1 

2 

1 1 
1 2 



1 

1 0.262433 

2 0.012774 

3 0.005155 
1 0.023223 
3 0.005575 

0.016024 
0.108998 
0.097751 
0.005624 
0.008794 
0.048875 
0.005624 
0.007209 
0.017229 
0.058045 
0.038323 
0.002354 
1 0.069670 
3 0.033448 
0.024036 
0.041514 
0.011247 
0.001433 
0.005624 
0.048875 
0.022494 
0.014417 
0.032937 
0.003923 
0.025773 
0.027873 
0.040060 



Assume spatial variations with wavelengths 2ir/k that 
are much larger than the effective collision range R e s- In 
this limit, all fields including the distribution function / 
and the density p are constant inside a circle which is 
centered around position x and has radius R c s- The ex- 
ponential prefactor in Eq. (|18[) becomes small for radii 
Rj > R c ff and thus very effectively suppresses errors 
when / and p are crudely approximated far away from 
x. This allows us to approximate the value of the den- 
sity in the integral of Eq. (|14[) by p(x) for any radius 
Rj . Then the integrand is constant and the simple result 
Mj = 7ri?| +1J o(x) is obtained. Similarly, we can formally 
replace the distribution functions /(#i,Xi) by their value 
at the point x, /(0j,x), where i = 2,3, ••■ , M, in the 
collision integral (| 18|) . This effectively ignores field vari- 
ations within typical collision distances. After integrat- 
ing over the positions of all the collision partners, the 
equation becomes, 

f(6,x + TV,t +r) 

1 f n/2 d£ f - - ~ - 

-i / dSidS-, . . . dS M S{9 - f - *i) 
-n/2 V J 

(101) 



p(x) A/ -! 

/(fl 1) x)/(fl 2 ,x).../(fl M) x). 



The approximative collision term on the r.h.s. contains 
only information from the point x but not from surround- 
ing points anymore. In an Enskog equation, it is the dif- 
ferences in the field values around point x, which account 
for the so-called collisional momentum transfer. There- 
fore, we have effectively removed this transfer, and Eq. 
(|101l) can be seen as the Boltzmann limit of the more 
general Enskog-like kinetic equation, Eq. (|T5)) . More- 
over, for M = 2, Eq. (| 10 1[) can be directly compared 
with the equation postulated by Peshkov et al. [21|. By 
rewriting Eq. ([lOljl into the format of Eq. (|47|) in or- 
der to investigate the linear stability, one sees that the 
original integrand F(0\, 02, ■ ■ ■ &m) is replaced by the fol- 
lowing, 



(102) 



M 



1 — (M - 1) 



S90 
9o 



i=2 



TABLE II: Integrals for M = 3 denned in Appendix C. 



Again, nonlinear perturbations of order 0{5 2 ) were 
dropped in this expansion. Comparing this with Eq. 
(|49|) . we note that it can be formally obtained by set- 
ting z = in the original kernel. Using the dimension- 
less wavenumber k = Xk and the definition of effective 
collision radius, Eq. (J3|), we have z = —R 2 s k 2 /(4M\ 2 ). 
Thus, for z to approach zero, either A = X/R c g must 
go to infinity or the wave number must be zero. This 
confirms that the approximative Boltzmann-likc equa- 
(I101D and hence the kinetic equation proposed in 



tion 
Rcf. 



present the zero wave number or infinite A 
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A=8 
- A=l 
■ A=l/4 



I I I I I I 3aJI I I I U 1_1 U 

1 2 3 1 2 3 

FIG. 17: Real part of the growth rate w as a function of 
wave number for M — 2, at r/ = 0.99t?c and with density 
po = 0.1. Part a) shows the first two longitudinal modes, 
whereas the first transversal mode is shown in b). The re- 
sults of the Enskog-like equation, Eq. (|47|) . for various ratios 
of mean free path to effective radius, A, are shown in dot- 
ted, dot-dashed, and dashed lines, whereas the results for the 
Boltzmann approach, Eq. (|101|) are given by gray solid lines. 
For the case where A is as large as 8, the curves of the origi- 
nal Enskog-like equation collapse with those of the Boltzmann 
approximation. 



limit where collisional momentum transfer is irrelevant. 
For M = 2 we compare the stability of the Enskog and 
the approximative Boltzmann equations in Fig. [T7I where 
the growth rates for perturbations of the ordered state are 
plotted. As expected we observe that both approaches 
agree exactly for zero k and very large A > 8. When 
A is decreased, collisional momentum transfer becomes 
more important and the growth rates become more neg- 
ative compared to the transfer- free case. Especially for 
A = 1/4 the difference is very pronounced, even at small 
Xk < 1. Hence, collisional momentum transfer makes the 
ordered phase even more stable. This is consistent with 
results for other systems [2(|, where collisional momen- 
tum transfer leads to an additional contribution to the 
viscosity and thermal conductivity, causing a stronger at- 
tenuation of sound and shear modes. In Fig. Q~7]one also 
sees that the sound mode is not as strongly affected by 
collisional momentum transfer as the other modes. 
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